#!/usr/bin/env python

###############################################################################
#
#
###############################################################################

#
# MODULES
#
import sys
import os
import re
import glob
import string

#from numarray import *
from pylab import *
#import Numeric

#
# ARGUMENTS
#

dir = 'monthly'

if len(sys.argv) < 2:
    print 'Usage: %s <dir>' % sys.argv[0]
    raise SystemExit
else:
   dir = sys.argv[1]


thresh_maxwind = 45.0
thresh_minpsl = 980.0

runs = [2, 3, 4, 5, 6, 7, 8, 10, 11, 13, 14, 15, 17, 18]
months = [7, 8, 9, 10, 11, 12]

past_tccounts = dict()
future_tccounts = dict()

pf = 'past'
for month in months:
    past_tccounts[month] = []
    for nrun in runs:
        count = 0
        fname = '%s/traj_%d_%s_summary_%d' % (dir, nrun, pf, month)
        if os.path.exists(fname):
            for line in open(fname).readlines():
                # test wind/psl criteria
                elems = string.split(line)
                if not elems[0].startswith('ID'): # ignore header
                    if (float(elems[7]) >= thresh_maxwind) & (float(elems[8]) <= thresh_minpsl):
                        count = count + 1
        past_tccounts[month].append(count)


pf = 'future'
for month in months:
    future_tccounts[month] = []
    for nrun in runs:
        count = 0
        fname = '%s/traj_%d_%s_summary_%d' % (dir, nrun, pf, month)
        if os.path.exists(fname):
            for line in open(fname).readlines():
                # test wind/psl criteria
                elems = string.split(line)
                if not elems[0].startswith('ID'): # ignore header
                    if (float(elems[7]) >= thresh_maxwind) & (float(elems[8]) <= thresh_minpsl):
                        count = count + 1
        future_tccounts[month].append(count)

#print past_tccounts
#print future_tccounts

# compute means and standard errors
past_means = dict()
future_means = dict()


# plot
figure(1)
hold(True)

for month in months:
    past_means = mean(past_tccounts[month])
    future_means = mean(future_tccounts[month])

    width = .4
    bar(month - width*.5, past_means, width*.5,0,'r')
    bar(month, future_means, width*.5,0,'b')

    legend(('Past','Future'))
    t = 'Monthly Storm Counts (Max Wind >= %2.2f; Min PSL <= %2.2f)' % (thresh_maxwind, thresh_minpsl)
    title(t)

    xlabel('Month')
    ylabel('Storm Count')

    xticks(months)
    
show()
